Adaptive molecular resolution via a continuous change of the phase space 
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For the study of complex synthetic and biological molecular systems by computer simulations 
one is still restricted to simple model systems or to by far too small time scales. To overcome this 
problem multiscale techniques are being developed for many applications. However in almost all 
cases, the regions of different resolution are fixed and not in a true equilibrium with each other. We 
here give the theoretical framework for an efficient and flexible coupling of the different regimes. 
The approach leads to an analog of a geometry induced phase transition and a counterpart of the 
equipartition theorem for fractional degrees of freedom. This provides a rather general formal basis 
for advanced computer simulation methods applying different levels of resolution. 

PACS numbers: 05.10.-a, 05.20.-y, 02.70.-c 



A long standing and often most challenging problem 
in condensed matter physics, up to some simple crys- 
talline materials, is to understand the microscopic origin 
of macroscopic properties. While in certain cases a micro- 
scopic scale can be clearly separated from a macroscopic 
one, this is not the case for most experimental systems 
- details of the local interaction and generic/universal 
aspects are closely related. Already the proper determi- 
nation of the fracture energy and crack propagation in 
crystalline materials requires a hierarchical and interre- 
lated description, which links the breaking of the inter- 
atomic bonding in the fracture region to the response of 
the rest of the system on a micron scale The pres- 
ence of microscopic chemical impurities in many metals 
and alloys changes their macroscopic mechanical behav- 
ior Even more complicated are synthetic and bio- 
logical soft matter systems. Whether one is dealing with 
the morphology, the glass transition of a polymeric sys- 
tem, the function of a molecular assembly, e.g. for elec- 
tronic applications or studies ligand-protein recognition 
or protein-protein interaction, in all cases the generic soft 
matter properties, such as matrix or chain conformation 
fluctuations, and details of the local chemistry apply to 
roughly the same length scales. For all these problems, 
which due to their complexity are heavily studied by com- 
puter simulations, there is a common underlying physical 
scenario: the number of degrees of freedom (DOFs) in- 
volved is very large and the exhaustive exploration of the 
related phase space is prohibitive. For many questions, 
however, such a deep level of detail in the description is 
only required locally. 

Theoretical methods employed to study these sys- 
tems span from quantum-mechanical to macroscopic 
statistical approaches. Their efficiency and scope in- 
creases significantly if two or more such different ap- 
proaches are combined into hybrid multiscale schemes. 
This is the case for the quantum based QM/MM ap- 
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proach M and that of dual scale resolution techniques 

EillTiillSllllIllIllIllIllIllIllIl aim- 
ing at bridging the atomistic and mesoscale length scale. 
However, the common feature and limitation of all these 
methods is the fact that the regions or parts of the system 
treated at different level of resolution are fixed and do 
not allow for free exchange. It is exactly this constraint 
one has to overcome in order to study typical complex, 
fluctuating molecular systems with a higher computa- 
tional efficiency. What is needed is an approach that 
allows to zoom into a specific area which, even though 
the number of DOFs treated is different, stays in equi- 
librium with the more coarse grained surrounding. The 
requirement above has been fulfilled by the recently in- 
troduced AdResS scheme [19|, |2Qj, but is more general 
as the above examples show. In computational terms 
such a simulation translates into a scheme of changing 
the number of DOFs on demand and on the fly in a se- 
lected region. While examples of robust numerical algo- 
rithms have been given p!9l l20l |21| the underlying concept 
still requires to be put into a solid, rigorous theoreti- 
cal framework, making it applicable also for the other 
schemes mentioned before. This is the subject of the 
present Letter. The continuous transition from a less 
to a more coarse grained description (and vice-versa) is 
explained in terms of a "geometrically induced first or- 
der phase transition", where the similarity with a stan- 
dard phase transition is further put into the context of 
non-integer dimensions of the phase space as DOFs are 
slowly switched off/on. In this way the present theory 
leads to the result of generalizing the equipartition theo- 
rem to non-integer dimensions and shows how to obtain 
relevant thermodynamic quantities within a continuous 
variable resolution of phase space. In this sense, the orig- 
inal problem of high dimensional systems, multiscale in 
nature, can be considerably simplified on the basis of 
a general and rigorous statistical mechanics framework. 
Because of the generality of the theoretical argument, 
this should be applicable to other multiscale simulation 
approaches. 

Without loss of any generality let us assume a system 
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of molecules in a volume V, modeled on a rather coarse- 
grained level. Now let us further assume that in a certain 
subvolume V a higher resolution is needed, i.e., to study 
some function. This is a typical situation, which one 
encounters in many systems, e.g., proteins or functional 
molecular assemblies. In statistical terms this translates 
into saying that the resolution employed in one region is 
lower (or higher) than in the rest. Thus the number of 
DOFs of the molecule in such a region is lower (or higher). 
For simplicity, we divide the volume V into two parts A 
and B. In region A, each molecule is characterized by ua 
DOFs, and in region B by ub- For example, in region 
A one has higher resolution and a molecule can be con- 
sidered as a collection of atoms linked by springs while 
in B the molecule consists only of its center of mass and 
a spherical excluded volume as in our previous numeri- 
cal example |l9l |2C| . The natural question now is how to 
reach true thermodynamic equilibrium between the two 
regions with the same overall structure of the system on 
both sides? At t his p oint we assume, tested numeric ally 
in previous work[ll l^illSlllllllllllllilEg 
that for a state point {p,T) it is in general possible to 
reduce the many body potential of the higher resolution 
representation into a dimensionally reduced effective po- 
tential. The latter, when applied to a system composed 
exclusively of molecules with the lower resolution, repro- 
duces the statistical properties of a system composed ex- 
clusively of highly resolved molecules, when analyzed ac- 
cordingly. However, treating one overall system with re- 
gions of different resolution, requires a special attention. 
The problem of the changing DOFs must be addressed 
in a way that A and B are in thermodynamical equilib- 
rium with each other and, additionally, the general over- 
all structure is the same. To address this we visualize 
the free energy F of the system as a function of the posi- 
tion X. (In order to do so we can divide the system into 
" large enough" equal slabs) . The free energy is a thermo- 
dynamic potential and hence in thermodynamic equilib- 
rium F{x) = Fa, constant in region A, and F{x) = Fb, 
constant in region B. Fa and Fb are in general different 
since F is an extensive quantity and ua ^ n n (29iJ . Note, 
however, that despite this free energy difference, which 
stems exclusively from the different levels of a molecular 
representation in regions A and B, the chemical poten- 
tials must be equal in both regions. This guarantees that 
the molecules experience no spurious driving force which 
would pull them from one region into the other due to 
the choice of the level of resolution. Figure sketches a 
"typical" free energy profile across the system. We focus 
now on the transition regime A between the two subsys- 
tems, i.e., —d < X < +d, where the points —d and +d 
denote the boundaries of A with the regions A and B, 
respectively, see figure H In A we gradually change the 
level of resolution and consequently the value of F{x). 
The width of A is set by the range of the effective pair 
potential between molecules Ts^- Our system is in equi- 
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FIG. 1: The free energy F as a function of a;. A and B are 
the regions with high and low levels of detail, respectively, 
while A is a transition regime. The constant values (minima) 
of Fa and Fb are arbitrary and do not play any role in the 
following treatment. 



librium, which implies that at the boundaries 
dFA{x) 
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If this condition did not hold, a molecule would 'see' a 
free energy gradient along x within the same level of 
resolution leading to a drift along the x axis. Next, 
1^1^^ and the same for B, 

Here Na and Nb are the numbers 



let us write 

OFb _ OFb BNb Que 
dx dNs duB dx 

of molecules in A and B, respectively. ON a /Qua and 
ONb I dns are two non-zero constants, while dFA/dNA — 
fiA and OFb/ONb = y^B, where ^ is the chemical poten- 
tial. Note that due to equilibrium /i^ = A's 7^ 0. The 
condition of Eq^Jis hence reduced to: 
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= 0. 



(2) 



Thus, we can formally describe the switching on and off 
of a given DOF via a weighting function w{x) such that 
w{x) = l;'ix G A and w{x) = 0;Vx G B, with zero slope 
at the boundaries of A. In accordance with Eq[21 this 
requires w{x) to be continuous up to the first derivative 
and monotonically goes from the value one to zero in 
the region A. The latter requirement reflects the fact 
that we would like to switch gradually on-the-fly with- 
out extra equilibration from more DOFs to less (or vice- 
versa) js^- The switching procedure implies that in the 
transition regime, where < w{x) < 1, we deal with 
fractional DOFs, i.e., by switching on/off a DOF we con- 
tinously change the dimensionality of the phase space. To 
rigorously describe the fractional phase space we resort 
here to fractional calculuSjSl, _33, _34, J5j. Accord- 
ing to Refs.js^, 01 I35II and to the formula for dimen- 
sional regularization '36'|, the infinitesimal volume ele- 
ment of the fractional configurational space is defined 
as dVa = d"a;r(a/2)/(W2r(a)) = \x\'^-^dx/r{a) = 
dx" / (aT{a)) where the positive real parameter a de- 
notes the order of the fractional coordinate differential 
and F is the gamma function. To make the connection 
with the switching on/off of DOFs we consider in our 
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case each a as a value that ■w{x) can take and apply the 
formalism to each DOF separately. Hence, we have to 
go beyond Refs. 0, 1^3 so that each DOF can take its 
own value of parameter a according to the level to which 
the particular DOF is switched on. Such a formulation 
leads to an interpretation of the switching as a " geomet- 
rically induced phase transition" . Here, we deal with a 
representation-driven geometrical transition, which has 
formal similarities with the standard concept of a first- 
order phase transition|37|. In fact, while in a standard 
first-order phase transition we have a latent heat, e.g. of 
solidification, in this case we have a latent heat due to the 
fact that the free energy density becomes position depen- 
dent and we need to furnish or remove " latent" heat from 
an external bath to compensate the free energy gradient 
due to changing the number of DOFs in A. 

To explain this concept in more detail we define the 
temperature in A in a usual way as Ta — 2 {KA)/nA, 
in B, Tg = 2 (Kb) /riB, where {Ka) and {Kb) are the 
average kinetic energies of a molecule in regions A and 
B, respectively; since we have equilibrium, Ta = Tb- In 
a similar way we would also like to define the tempera- 
ture in the interface region A as Ta = 2 {KA)/nA with 
(Ka) and being the average kinetic energy and the 
number of DOFs of a molecule in a given slab x in A, 
respectively. However, for such a definition we first need 
to determine how and riA scale with w{x). Let us 
demonstrate this with a simple example. In figure |21 we 
have a disc representing a 2-d molecule with three DOFs 
in the high resolution region A, i.e., two translational 
DOFs of the center of mass R = [Rx, Ry) and one rota- 
tional DOF around the center of mass characterized by 
angle O, the transition representation of the molecule in 
the region A, and the coarse grained molecule in the re- 
gion B with only two translational DOFs. The kinetic 




FIG. 2: The molecular resolution of a simple 2-d circular 
molecule in the high and low resolution regions A and B and 
the transition region A, respectively. 

energy in the region A, setting the mass, the molecule's 
radius, and the Boltzmann constant m = r = kB ~ I, 
is: Ka = [R^ -l-p|]/2 = [R.J^ + Ry^ + e'^]/2, and in re- 
gion B: Kb = R^/2 = [r/ + Ry\l2. Here pe denotes 
the angular momentum. According to the equipartition 
theorem each full quadratic DOF, i.e., R,J,^ Ry, and O, 
contributes to the kinetic energy with an amount of T/2. 



Hence, Ta = 2 < Ka >/3, Tb = 2 < Kb >/2. In the re- 
gion A however, Q is not a full DOF, because its weight 
varies between zero and one. This has to be considered, 
when calculating the local temperature. To determine 
such a contribution one should account for the follow- 
ing: Let us consider a given value w{x) — a and use 
a as a variable parameter. The instantaneous kinetic 
energy of a fractional DOF should gradually vanish as 
the DOF is slowly switched off (or vice- versa), thus in 
the region A, our ansatz for the kinetic energy associ- 
ated to 9 is f{a)pQ/2, where /(a) is a monotonic func- 
tion in A with /(I) = l;/(0) = 0. Apart from these 
requirements we do not need to specify the exact form 
of /(a) at this point. Accordingly, the kinetic energy 
in the position with the coordinate x in the region A 

• 2 ■ 2 

is Ka — [Rx + Ry + f{a)pQ]/2. For the fractional 
quadratic DOF O we can then write the partition func- 
tion as: 

exp(-/3F„) J exp{-pf{a)pl/2)dVa = 



2C / exp(-/3/(a)p4/2) \pe 
Jo 

2"/2CT(a/2 



r(a) 



r{a) 



(3) 



where C is a normalization constant, (3 — 1/T, and Fa 
the free energy associated with the fractional DOF O, 
respectively |38j . The consequence of Eq|31 is the frac- 
tional analog of the equipartition theorem: 



djPFq 
d/3 



a 



aT 



(4) 



where (Ka) is the average kinetic energy per fractional 
quadratic DOF with the weight a. Thus, for a = 0, 1 
we obtain the correct limits in the coarse-grained and 
fully resolved regimes, respectively, with the correct con- 
tributions to the kinetic energy. Furthermore, we have 
Ta = 2 (Ka) /ua = (2 -t- a)T/nA- To satisfy the equi- 
librium condition: Ta = Tb = Ta = T we must set 
TT-A — 2 + a, which is in accordance with the "intuitive" 
definition in Ref. 0|. The number of quadratic DOFs 
and the average kinetic energy thus scale linearly with 
w{x). Eq0|also tells us that, although the equipartition 
is independent of the specific choice of f{a), since the av- 
erage kinetic energy scales as a also the instantaneous one 
should scale in the same way. This means that /(a) = a 
and it is determined by the fractional character of the 
phase space. Note that for non-quadratic DOFs the func- 
tional form of /(a) is generally more complicated [l^, 
however not needed for the present purpose. 

Based on this theoretical considerations we have re- 
cently derived an efficient particle-based MD simulation 
scheme ^3 For intermolecular force calculation we 
use an interpolation formula for the force acting between 
centers of mass of given molecules a and /3: 

= w{xa)w{xp)¥'^J^-^ + [1 - w{xa)w{xp)]F''^f, (5) 
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where Xa and Xji are the center of mass coordinates of the 
molecules a and /?, respectively, FJ^*^™ is the sum of all 
pair atom interactions between explicit atoms of molecule 
a and explicit atoms of molecule (3 and F^^ is the total 
force between the centers of mass of the respective two 
molecules l2(ij. This ansatz satisfies Newton's Third 
Law and takes into account the transfer of the turned off 
explicit DOF onto the molecular center of mass,40j. Each 
time a given molecule crosses a boundary between differ- 
ent regimes it gains or looses (depending on whether it 
leaves or enters the region B) its equilibrated rotational 
DOF while retaining its linear momentum. By extension 
of the equipartition theorem to fractional DOFs we are 
able to define the means to supply the latent heat, which 
is required or removed for the transition from coarse- 
grained to high resolution or vice-versa. Since this gen- 
eralized equipartition theorem also applies to the frac- 
tional quadratic DOFs standard thermostats based on 
the fluctuation-dissipation theorem are applicable. 



In conclusion, we provided the statistical mechanics 
foundation for an efficient computational scheme that 
concurrently couples different length scales via different 
levels of resolution, i.e., atomistic and mesoscopic length 
scales, by adapting the resolution on demand. The tran- 
sition region is well defined by the here introduced gen- 
eralization of the equipartition theorem for fractional di- 
mension of phase space. While it directly applies to 
a scheme recently tested by the authors [1^ it in 
the same way should also provide the general theoreti- 
cal framework to extend other commonly used schemes, 
such as [T^ ITgL It^ towards a truly adaptive multiscale 
simulation scheme. 

We thank M. Deserno, D. Andrienko, and L. Ghir- 
inghelli for helpful comments on the manuscript. This 
work is supported in part by the Volkswagen founda- 
tion. M. P. acknowledges the support of the Ministry 
of Higher Education, Science and Technology of Slovenia 
under grant No. Pl-0002. 



[1] J. Mottler, S. Barsky, and M. O. Robbins, Phys. Rev. 

Lett. 89, 148304 (2002). 
[2] D. E. Jiang and E. A. Carter, Acta Materialia 52, 4801 

(2004) . 

[3] G. Lu and E. Kaxiras, Pliys. Rev. Lett. 94, 155501 

(2005) . 

[4] A. Laio, J. VandeVondele, and U. Rothlisberger, J. 

Chem. Pliys. 116, 6941 (2002). 
[5] L. Delle Site, C. F. Abrams, A. Alavi, and K. Kremer, 

Pliys. Rev. Lett. 89, 156103 (2002). 
[6] L. Delle Site, S. Leon, and K. Kremer, J. Am. Chem. 

Soc. 126, 2944 (2004). 
[7] E. Villa, A. Balaeff, L. Mahadevan, and K. Schulten, 

Multiscale Model. Simul. 2, 527 (2004). 
[8] M. Neri, C. Anselmi, M. Cascella, A. Maritan, and 

P. Carloni, Phys. Rev. Lett. 95, 218102 (2005). 
[9] H. Rafii-Tabar, L. Hua, and M. Cross, J. Phys.: Condens. 

Matter 10, 2375 (1998). 
[10] J. Q. Broughton, F. F. Abraham, N. B. Bernstein, and 

E. Kaxiras, Phys. Rev. B 60, 2391 (1999). 
[11] J. A. Smirnova, L. V. Zhigilei, and B. J. Garrison, Comp. 

Phys. Commun. 118, 11 (1999). 
[12] G. Csanyi, T. Albaret, M. C. Payne, and A. D. Vita, 

Phys. Rev. Lett. 93, 175503 (2004). 
[13] S. T. O'Connell and P. A. Thompson, Phys. Rev. E 52, 

R5792 (1995). 

[14] N. G. Hadjicostantino, Phys. Rev. E 59, 2475 (1999). 
[15] J. Li, D. Liao, and S. Yip, Phys. Rev. E 57, 7259 (1998). 
[16] E. G. Flekkoy, G. Wagner, and J. Feder, Europhys. Lett. 
52, 271 (2000). 

[17] R. Delgado-Buscalioni and P. V. Coveney, Phys. Rev. E 

67, 046704 (2003). 
[18] J. Baschnagel et al.. Adv. Polym. Sci. 152, 41 (2000). 
[19] M. Praprotnik, L. Delle Site, and K. Kremer, J. Chem. 

Phys. 123, 224106 (2005). 
[20] M. Praprotnik, L. Delle Site, and K. Kremer, Phys. Rev. 

E 73, 066701 (2006). 
[21] C. F. Abrams, J. Chem. Phys. 123, 234101 (2005). 
[22] R. L. Henderson, Phys. Lett. 49A, 197 (1974). 



[23] A. P. Lyubartsev and A. Laaksoneri, Phys. Rev. E 52, 

3730 (1995). 
[24] A. K. Soper, Chem. Phys. 202, 295 (1996). 
[25] W. Tschop, K. Kremer, O. Hahn, J. Batoulis, and 

T. Biirger, Acta Polym. 49, 75 (1998). 
[26] D. Reith, M. Piitz, and F. Miiller-Plathe, J. Comput. 

Chem. 24, 1624 (2003). 
[27] S. Izvekov, M. Parrinello, C. B. Burnham, and G. A. 

Voth, J. Chem. Phys. 120, 10896 (2004). 
[28] S. Izvekov and G. A. Voth, J. Chem. Phys. 123, 134105 

(2005). 

[29] The fact that Fa and Fb need not to be equal illustrates 
the limitations of the analogy to phase transitions, where 
the free energy at the transition point is the same in the 
two phases. 

[30] Eq.2 expresses the fact that we cannot change ua and ub 
otherwise than by changing Na and Nb, respectively. 
In contrast, we also change the number of DOFs in A 
by changing the number of DOF per molecule with the 
weighting function w^x). 

[31] T. F. Nonnenmacher, J. Phys. A: Math. Gen. 23, L697S 
(1990). 

[32] R. Hilfer, editor. Applications of fractional calculus in 
physics, World Scientific Publishing, Co. Pte. Ltd., Sin- 
gapore, 2000. 

[33] K. Cotrill-Shepherd and M. Naber, J. Math. Phys. 42, 
2203 (2001). 

[34] V. E. Tarasov, Chaos 14, 123 (2004). 

[35] V. E. Tarasov, Phys. Rev. E 71, 011102 (2005). 

[36] J. C. Collins, Renormalization, Cambridge University 
Press, Cambridge, 1984. 

[37] R. M. Lyndell-Bell, Mol. Phys. 86, 1353 (1995). 

[38] R. P. Feynman and A. R. Hibbs, Quantum mechanics 
and path integrals, McGraw-Hill, Inc., New York, 1965. 

[39] M. Deserno, private communication. 

[40] One could, in principle, also define the coupling scheme 
using the respective potentials instead of forces as f/a/3 = 
w{xc)w{xp)U^f"' + [l-w{xc.)w{x^)]U^'lj. Note that such 
a linear combination depends on the absolute positions 



5 

of the particles and not only on their relative distances. obey Newton's Third Law unless w{x) = const. 

A force defined as a gradient of this potential would not 



